//
//  Obstacle.cpp
//  NSESolver
//
//  Created by Patrick Zulian on 5/5/11.
//  Copyright 2011 ICS - USI - Lugano. All rights reserved.
//


#include "Obstacle.h"
#include "CFDProblem.h"


namespace nsesolver {
    
    
    
    Obstacle::Obstacle(const int imax, const int jmax)
    : cells(), obsCellMap(imax + 2, jmax + 2)
    { 
        obsCellMap.allSet(false);
        obsCellMap(Range(0), Range(0, jmax+2)) = true;
        obsCellMap(Range(imax+1), Range(0, jmax+2)) = true;
        obsCellMap(Range(0, imax+2), Range(0)) = true;
        obsCellMap(Range(0, imax+2), Range(jmax+1)) = true;
       
    }
    
    Obstacle::~Obstacle()
    { }
    
    
    void 
    Obstacle::setBoundary(CFDProblem &problem)
    {
        Matrixd & u(problem.U());
        Matrixd & v(problem.V());
        
        std::vector< obstacle_cell >::iterator it;
        for (it = cells.begin(); it != cells.end(); ++it) {
            int i = it->i, 
            j = it->j;
            switch (it->type) {
                    
                case ONorth: { 
                    u(i-1, j) = -u(i-1, j+1);
                    u(i, j) = -u(i, j+1);
                    
                    v(i, j) = 0;
                    break;
                }
                    
                case OSouth: { 
                    u(i-1, j) = -u(i-1, j-1);
                    u(i, j) = -u(i, j-1);
                    
                    v(i, j-1) = 0;
                    break;
                }
                    
                case OWest: { 
                    u(i-1, j) = 0;
                    
                    v(i, j-1) = -v(i-1, j-1);
                    v(i, j) = -v(i-1, j);
                    break;
                }
                    
                case OEast: { 
                    u(i, j) = 0;
                    
                    v(i, j-1) = -v(i+1, j-1);
                    v(i, j) = -v(i+1, j);
                    break;
                }
                    
                case ONorthWest: { 
                    u(i-1, j) = 0;
                    u(i, j) = -u(i, j+1);
                    
                    v(i,j) = 0;
                    v(i, j-1) = -v(i-1, j-1);
                    break;
                }
                    
                case OSouthEast: { 
                    u(i, j) = 0;
                    u(i-1, j) = -u(i-1, j-1);
                    
                    v(i, j-1) = 0;
                    v(i, j) = -v(i+1, j);
                    break;
                }
                    
                case ONorthEast: { 
                    u(i, j) = 0;
                    u(i-1, j) = -u(i-1, j+1);
                    
                    v(i, j) = 0;
                    v(i, j-1) = -v(i+1, j-1);
                    break;
                }
                    
                case OSouthWest: { 
                    u(i-1, j) = 0;
                    u(i, j) = -u(i, j-1);
                    
                    v(i, j-1) = 0;
                    v(i, j) = -v(i-1, j); 
                    break;
                }
                    
                case ONull: { 
                    u(i, j) = v(i, j) = 0; //TODO
                    u(i-1, j) = v(i, j-1) = 0; //TODO
                    break;
                }
                default:
                    break;
                    
            }
            
        }
    }
    
    void 
    Obstacle::setPBoundary(CFDProblem &problem)
    {
        Matrixd & p(problem.P());
        const double dx2 =problem.DX() * problem.DX(),
        dy2 = problem.DY() * problem.DY();

        std::vector< obstacle_cell >::iterator it;
        for (it = cells.begin(); it != cells.end(); ++it) {
            int i = it->i, 
            j = it->j;
            
            switch (it->type) {
                    
                case ONorth: { 
                    p(i, j) = p(i, j+1);
                    break;
                }
                    
                case OSouth: { 
                    p(i, j) = p(i, j-1);
                    break;
                }
                    
                case OWest: { 
                    p(i, j) = p(i-1, j);
                    break;
                }
                    
                case OEast: { 
                    p(i, j) = p(i+1, j);
                    break;
                }
                    
                case ONorthWest: { 
                    p(i, j) = ( dx2 * p(i, j+1) + dy2 * p(i-1, j) )/(dx2 + dy2);
                    break;
                }
                    
                case OSouthEast: { 
                    p(i, j) = ( dx2 * p(i, j-1) + dy2 * p(i+1, j) )/(dx2 + dy2);
                    break;
                }
                    
                case ONorthEast: { 
                    p(i, j) = ( dx2 * p(i, j+1) + dy2 * p(i+1, j) )/(dx2 + dy2);
                    break;
                }
                    
                case OSouthWest: { 
                    p(i, j) = ( dx2 * p(i, j-1) + dy2 * p(i-1, j) )/(dx2 + dy2); 
                    break;
                }
                    
                case ONull: { 
                    p(i, j) = 0; //TODO
                    break;
                }
                default:
                    break;
                    
            }
            
        }
    }
    
    
    void 
    Obstacle::setFGBoundary(Matrixd & F, Matrixd & G, const CFDProblem &problem)
    {
        const Matrixd & u(problem.U());
        const Matrixd & v(problem.V());
        
        std::vector< obstacle_cell >::iterator it;
        for (it = cells.begin(); it != cells.end(); ++it) {
            int i = it->i, 
            j = it->j;
         
            switch (it->type) {
                    
                case ONorth: { 
                    G(i, j) = v(i, j);
                    break;
                }
                    
                case OSouth: { 
                    G(i, j-1) = v(i, j-1);
                    break;
                }
                    
                case OWest: { 
                    F(i-1, j) = u(i-1, j);
                    break;
                }
                    
                case OEast: { 
                    F(i, j) = u(i, j);
                    break;
                }
                    
                case ONorthWest: { 
                    F(i-1, j) = u(i-1, j);
                    G(i, j) = v(i, j);
                    break;
                }
                    
                case OSouthEast: { 
                     G(i, j-1) = v(i, j-1);
                     F(i, j) = u(i, j);
                    break;
                }
                    
                case ONorthEast: { 
                    F(i, j) = u(i, j);
                    G(i, j) = v(i, j);
                    break;
                }
                    
                case OSouthWest: { 
                     G(i, j-1) = v(i, j-1);
                     F(i-1, j) = u(i-1, j);
                    break;
                }
                    
                case ONull: { 
                    G(i, j) = u(i, j);
                    F(i, j) = v(i, j);
                    break;
                }
                default:
                    break;
                    
            }
            
        }
    }
    
    void
    Obstacle::assignFlags()
    {
        std::vector< obstacle_cell >::iterator it;
        for (it = cells.begin(); it != cells.end(); ++it)
        {
            int neigh=0;
            bool west=false, east=false, south=false, north=false;
            const int i=it->i, j=it->j;
            
            if(isObstacleCell(i-1, j))
                ++neigh;
            else
                west=true;
            
            if(isObstacleCell(i+1, j))
                ++neigh;
            else
                east=true;
            
            if(isObstacleCell(i, j-1))
                ++neigh;
            else
                south=true;
            
            if(isObstacleCell(i, j+1))
                ++neigh;
            else
                north=true;
            
            if(neigh==1){
                NSE_SOLVER_ASSERT(false);
                std::cout<<"Cell "<<i<<","<<j<<" has only one neighboring obstacle cell"<<std::endl;
                exit(0);
            }
            
            if (neigh==4) {
                it->type=ONull;
                continue;
            }
            
            if(neigh==3){
                if(west)
                    it->type=OWest;
                else if(east)
                    it->type=OEast;
                else if(south)
                    it->type=OSouth;
                else if(north)
                    it->type=ONorth;
                
                continue;
            }
            
            if(neigh==2){
                if(west && south)
                    it->type=OSouthWest;
                
                if(west && north)
                    it->type=ONorthWest;
                
                if(east && south)
                    it->type=OSouthEast;
                
                if(east && north)
                    it->type=ONorthEast;
                
                continue;
            }
        }
       
        //setContactBoundaryCells();
    }
    
    
    void 
    Obstacle::addCell(int i, int j)
    {
        obstacle_cell cell = { i+1, j+1, ONull };
        cells.push_back(cell);
        obsCellMap(i+1, j+1) = true;
    }
    
}
